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Abstract. The final burst of gravitational radiation emitted by coalescing binary 
neutron stars carries direct information about the neutron star fluid, and, in particular, 
about the equation of state of nuclear matter at extreme densities. The final merger 
may also be accompanied by a detectable electromagnetic signal, such as a gamma-ray 
burst. In this paper, we summarize the results of theoretical work done over the past 
decade that has led to a detailed understanding of this hydrodynamic merger process 
for two neutron stars, and we discuss the prospects for the detection and physical 
interpretation of the gravity wave signals by ground-based interferometers such as 
LIGO. We also present results from our latest post-Newtonian SPH calculations of 
binary neutron star coalescence, using up to 10 6 SPH particles to compute with higher 
spatial resolution than ever before the merger of an initially irrotational system. We 
discuss the detectability of our calculated gravity wave signals based on power spectra. 

INTRODUCTION 

Coalescing binary neutron stars (NS) are among the most promising sources of 
gravitational radiation that should be detectable by future generations of gravity 
wave detectors. LIGO, VIRGO, GEO, and TAMA may ultimately not only serve 
to test the predictions of the theory of general relativity (GR), but could also yield 
important information on the interior structure of neutron stars, which cannot be 
obtained directly in any other way. 

Compact binary orbits decay through energy losses to gravitational radiation. 
So long as the separation between the two NS is large, the binary inspiral is well 
described by a point-mass treatment, modified to take into account finite-size and 
relativistic effects, which act only as small corrections. At the end of the inspiral, 
however, the process is inherently hydrodynamic in nature. Large tidal interactions 
can drive the system into dynamical instability, at which point a quasi-equilibrium 
treatment of the binary breaks down, and a numerical treatment is required to 
accurately model the system. 

Essentially all recent calculations agree on the basic picture that emerges for the 
final coalescence (see [1] and [2] for a complete list of references). As the dynamical 
stability limit is approached, typically at separations of r = 3 — 4:R NS , depending 
on the choice of parameters, the NS undergo a rapid radial plunge and merge in 



no more than a few rotation periods, much more quickly than would be predicted 
by a point-mass formula. In many cases, especially for binaries assumed to be 
initially synchronized, mass shedding sets in immediately after the stars first make 
contact. Material with a high specific angular momentum located in the outer re- 
gions of each NS is shed through the outer Lagrange points of the system, forming 
spiral arms that encircle the merger remnant left in the center as the NS cores 
merge. Eventually, the spiral arms also merge into a nearly axisymmetric torus 
around the dense inner core. For a stiff equation of state (EOS), a core with a 
significant ellipsoidal (triaxial) deformation can be maintained, and the configura- 
tion keeps radiating gravity waves well after the merger is completed. Softer EOS 
cannot support such a configuration stably, and any remnant produced will relax 
on a dynamical timescale toward a spheroidal (axisymmetric) configuration, which 
produces a negligible amount of gravity waves. 

The previous statements assume that the merger remnant formed is stable against 
gravitational collapse to a black hole. Unfortunately, Newtonian calculations are 
incapable of demonstrating such an effect. Post-Newtonian (PN) simulations can 
produce configurations that are unstable against collapse, but they are inherently 
unreliable because in conditions of strong gravity the basic assumptions of the PN 
expansion break down. Early full GR calculations indicate that merger remnants 
may very well be stable against collapse, so long as the EOS is assumed to be 
stiff enough [3]. The mass of the remnant should be nearly twice the mass of a 
single NS, which is generally taken to be M^s ~ 1-4 — 1.5M , and thus well over 
the maximum mass for a single, nonrotating NS. However, the remnants formed in 
binary coalescence are very rapidly and differentially rotating, which can increase 
the maximum stable mass to a much larger value [4,5]. 

BINARY NS COALESCENCE CALCULATIONS 

Nakamura and collaborators [6,7] were the first group to perform 3-D hydrody- 
namic calculations of binary NS coalescence, using an Eulerian grid-based code. 
Rasio and Shapiro [8] used the Lagrangian SPH (Smoothed Particle Hydrodynam- 
ics) method to calculate gravitational wave forms from binary coalescence events. 
Calculations performed since, using both Eulerian methods [10-14] and SPH [15-18] 
have focused on several aspects of the problem, including the effects of different ini- 
tial spins, mass ratios, NS EOS, and NS masses. Some groups have incorporated 
treatments of the nuclear physics involved in the merger [12-14,16-18] in order to 
study coalescing NS binaries as possible gamma-ray burst sources, and as possible 
birthplaces for r-process elements. 

Much of the early work on coalescing NS binaries assumed Newtonian gravity 
for simplicity. Later studies added a treatment of the radiation reaction, which 
is responsible for driving the system towards coalescence, either by adding a fric- 
tional drag term to model point-mass inspiral [15-18], or by an exact PN treat- 
ment [12-14]. In essence, 2.5PN radiation reaction terms (which scale like 1/c 5 ) 



are added onto a Newtonian framework, but all lower-order non-dissipative terms 
are ignored. Unlike adding a frictional drag term which dissipates energy according 
to the point-mass prediction, the lowest-order treatment of the radiation reaction 
allows for its effects to be included throughout the entire calculation, including the 
period after the merger remnant has formed. Unfortunately, however, Newtonian 
gravity is known to be a poor description of the physical problem at hand. Even 
NS with stiff EOS generate strong gravitational fields. During the final moments 
before merger, the velocities found in the system also become relativistic. Thus, 
the hydrodynamics of the actual coalescence can only be calculated properly by 
taking into account GR effects. 

The Newtonian limit also fails to describe accurately the onset of dynamical 
instability. PN effets combine nonlinearly with finite-size fluid effects and this can 
dramatically increase the critical binary separation (and thus lower the frequency) 
at which dynamical instability sets in. Indeed, the quasi-equilibrium description 
applies so long as 



where E equi i(r) is the energy of an equilibrium binary configuration at a given 
separation r. Equilibrium sequences have been calculated for both synchronized 
and irrotational binaries in Newtonian gravity [19], PN gravity [20], and recently in 
full GR [21,22] (see also Baumgarte, in this volume). It is generally found that the 
energy of NS binary configurations reaches a minimum at some critical separation, 
defining the innermost stable circular orbit (ISCO). Relativistic terms move the 
ISCO to larger separations, and also reduce the slope of the energy curve just 
outside the ISCO. Thus, the assumption of quasi-equilibrium, which is used to set 
up the initial configuration of the binary system, breaks down at a much larger 
separation than a Newtonian calculation would predict. In addition, NS binaries 
will already have developed significant infall velocities as they pass through the 
ISCO calculated for systems in strict equilibrium. 

Several groups have been working on full GR calculations of binary NS mergers, 
but only preliminary results have been reported so far [23]. Proving to be partic- 
ularly difficult is the extraction of wave forms from the boundaries of large 3-D 
grids, since extending the grids into the true wave zone would be too expensive 
computationally. The middle ground between the Newtonian treatments and full 
GR lies in PN hydrodynamic calculations of binary mergers. The authors [1,24], 
as well as Ayal et al. [25], have constructed a PN SPH code, described below, for 
calculating binary mergers. While it too serves only as an approximation to the 
proper physics which must go into a realistic calculation, the results do provide 
insight into the relativistic effects that simple Newtonian intuition fails to handle 
correctly. Additionally, they should serve as valuable checks for future full GR 
calculations, lending confidence to wave form predictions, and indicating to some 
degree the difference between real relativistic effects and numerical instabilities. 
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POST-NEWTONIAN SPH 



Our Post-Newtonian calculations use a formalism adapted from that of Blanchet, 
Damour, and Schaefer [26]. It includes all first-order (1PN) terms of GR, as well as 
the lowest-order radiation reaction terms (2.5PN). The latter are important because 
they provide the energy dissipation mechanism which drives the binary system to- 
ward dynamical instability and coalescence. Calculating various PN quantities 
requires the solution of seven additional Poisson equations for 1PN terms and an 
additional Poisson equation for the radiation reaction, all of which have compact 
support and thus can be solved by similar methods as the Newtonian gravitational 
potential. All hydrodynamic quantities in the formalism have been converted from 
a grid-based Eulerian approach to a particle-based Lagrangian approach, where 
particles represent distributions of matter in space defined by a spherically sym- 
metric smoothing kernel (rather than point-like objects). Pressure forces and other 
interactions between particles are handled by summing over neighboring particles. 
Poisson equations are solved by translating particle-based source term quantities 
onto a grid, and solving via FFT-based convolution methods [24]. 

Unfortunately, the consistent use of physically realistic NS parameters is impos- 
sible in this PN formalism. The compactness of a NS is given by 



and is typically assumed to fall near GM/Rc 2 pa 0.15 for a reasonable choice of 
EOS. Since several of the coefficients of the 1PN terms can be quite large, espe- 
cially for stiff EOS, we find that in many cases the 1PN corrections are larger than 
the Newtonian terms. Thus, we adapt the formalism, reducing the magnitude of 
the 1PN corrections by a factor of three (in effect decreasing the mass of the NS 
to Mipn = 0.5M Q ), but treating all the radiation reaction terms, of significantly 
smaller magnitude, at full strength. In essence, we employ a different speed of light 
for 1PN and 2.5PN terms. By comparison with result including radiation reaction 
but no 1PN terms whatsoever, we believe we can extrapolate in a qualitative sense 
toward the proper physically realistic case. By comparison, self-consistent PN cal- 
culations performed for NS with artificially small masses [7,25] have the advantage 
that they can be directly compared to full GR calculations, for which there can 
be no separation of relativistic terms into separate orders. Unfortunately, these 
calculations also suffer from a drastic and completely artificial reduction of the ra- 
diation reaction effects, which scale like M 2 5 . This produces a significant delay in 
the onset of dynamical instability past the point where it would be encountered for 
a physical set of parameters, and can lead to a qualitatively incorrect description 
of the subsequent merger. 





IRROTATIONAL BINARY COALESCENCE 



Our most detailed calculation performed to date uses N = 500, 000 particles per 
NS, corresponding to the highest spatial resolution ever for a binary coalescence cal- 
culation. The spatial resolution (smoothing length) achieved in the central regions 
of the stars is h 0.03Rns- The calculation was performed using an irrotational 
initial condition. This is generally thought to be the most realistic case since the 
viscous tidal locking timescale for two NS is expected to be considerably longer than 
the inspiral timescale [27]. Corotating (tidally locked) systems are motionless in a 
frame corotating with the binary, which allows for the use of relaxation techniques 
that can give a very accurate equilibrium initial condition (thus nearly completely 
eliminating spurious oscillations around equilibrium during the early phases of the 
dynamical evolution). Instead, for our irrotational calculation, we model the initial 
density and velocity profile of the NS as tidally stretched ellipsoids, with parame- 
ters drawn from the PN equlibrium calculations of Lombardi, Rasio, and Shapiro 
[20]. 

Since all NS in relativistic binary systems are expected to have masses that lie 
within a very narrow range [28], our calculation uses equal-mass NS. As the NS 
EOS is still poorly constrained, we choose a simple T = 3 polytropic EOS, i.e., the 
pressure is given in terms of the rest-mass density by P — kp^. Stiff EOS, such as 
this one, are capable of maintaining a long-lived ellipsoidal deformation after the 
binary merger, with a gravity wave signal that persists on a timescale much longer 
than the merger timescale [9]. 

Particle plots showing the evolution of the equal-mass irrotational binary system 
are shown in Fig. 1. We see that immediately prior to merger a large tidal lag 
angle develops. The inner edge of each NS leads the axis connecting the centers 
of mass of the binary components, and the outer regions lag behind. This effect 
is seen even in Newtonian simulations, but is greatly enhanced by the addition 
of 1PN correction terms. When first contact is made, a long vortex sheet forms 
at the interface between the two stars. Unlike the case of synchronized binaries, 
we do not see significant mass shedding from the outer Lagrange points of the 
system. The rotational speed of particles on the outer half of each NS is reduced 
in the irrotational case with respect to the synchronized case, and such particles 
remain bound and form the outer regions of the eventual merger remnant. At 
t = 25 we do see some hint of mass shedding, but not via the mechanism described 
above. Particles which have travelled the length of the vortex sheet and retained 
significant velocities end up being shed from the leading edge of the vortex sheet. 
It is important to note, though, that the amount of mass shed is extremely small, 
much less than 1% of the total mass, and that the velocities of the particles ejected 
are not sufficient to escape the gravitational potential of the remnant. We thus 
expect them to form an extremely tenuous halo around the central core. At late 
times, we see the formation of a remnant containing essentially all the mass that 
was orignally present in the system. As we are using a stiff EOS, the ellipsoidal 
deformation of the remnant is relatively large, and persists for late times after the 




FIGURE 1. Final merger of two identical T = 3 polytropes with an irrotational initial condition. 
SPH particles are projected onto the equatorial plane of the binary. The orbital rotation is 
counterclockwise. Spatial coordinates are given in units of the NS radius R. Times are given 
in units of the dynamical timescale of the system, which here is to = 0.07ms = 1. The orbital 
rotation is in the counterclockwise direction. 



merger. 

Density contours and velocity profiles in the equatorial plane of the binary are 
shown in Fig. 2. Velocities are shown in a frame corotating with the material. We 
see that the initially counterstreaming surfaces of the NS produce a vortex sheet, 
which is Kelvin-Helmholtz unstable. Large vortices develop along the surface of 
contact, both in the center of the forming merger remnant and also at a separation 
which seems to be roughly consistent with the misalignment of the leading edges of 
each NS, or about r = 0.5R. As the merger proceeds, we see that the vortices re- 
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FIGURE 2. Density contours and velocities along the equatorial plane in the corotating frame 
of the binary, for the same times as in Fig. 1. 

main coherent from t = 20 — 25, mixing material which was orignally located along 
the inner parts of each NS. All the while, the cores of the respective NS continue 
to inspiral toward the center of the merger remnant, until by t = 30 they have 
formed a single core, the vortices having merged together. This produces a char- 
acteristic differentially rotating pattern, with the center of the remnant spinning 
approximately twice as fast as the outer regions. 

GRAVITY WAVE SIGNALS AND SPECTRA 

We calculate the gravity wave signal for our mergers in the quadrupole approx- 
imation. The gravity wave strain h seen by an observer located a distance d from 



r=3, q=1.0, Trrotational 




10 20 30 

t 



FIGURE 3. Gravity wave signals calculated for coalescences with the same initial parameters 
but different numerical resolutions. The solid, dashed, and dot-dashed lines correspond to runs 
with 10 6 , 10 5 , and 10 4 SPH particles, respectively. 

the center of mass of the system along the rotation axis is given for the two polar- 
izations by 

c 4 (dh + ) = Q xx - Q yy (3) 
c\dh x ) = 2Q xy (4) 

where Q, the second time derivative of the quadrupole moment tensor, is given in 
SPH terms by 

Qij = E m^vf + x?vf + xfv?) (5) 
b 

where the summation is taken over all particles in the calculation. 

In Fig. 3, we show the gravity wave signals in both polarizations for the irro- 
tational run described above, as well as for runs with N = 50, 000 particles and 
iV = 5, 000 particles per NS. It is immediately apparent that the lowest resolu- 
tion run shows significant discrepancies from the other two, which agree with each 
other quite well over the entire time history of the merger. This is a welcome result, 
given that the vortex sheet appearing at the contact surface is Kelvin-Helmholtz 



unstable on all wavelengths, including those much smaller than our numerical res- 
olution. Calculations performed at different resolutions do show subtle differences 
in the exact location and size of the vortices. It is important to note, however, 
that it is the outer regions of the star, at lower density, that supply material to 
the vortex sheet. The high density cores of the two NS inspiral during the entire 
process, and provide the dominant component of the quadrupole moment and thus 
the gravity wave signal. The path traced out by the NS cores depends sensitively 
on gravitational forces and properties of the fluid, such as the EOS, but proves to 
be remarkably insensitive to the details of the flow in the "turbulent" boundary 
region. The conclusion to be drawn is that numerical convergence for a given set 
of initial conditions and physical assumptions is possible without requiring exces- 
sive computational resources, even for this difficult problem involving small-scale 
instabilities. 

Because the gravity wave signals expected to be seen from binary NS mergers 
are extremely close to the sensitivity limits of ground-based interferometers, it is 
important to identify which features in the power spectrum of the gravitational ra- 
diation will yield the most information on the physical parameters of NS. Following 
the method of Zhuge et al. [15], we compute the gravity wave power spectrum per 
unit frequency interval as 

^ = ~(^)f 2 (\h(f)\ 2 + \h,(f)\ 2 )- (6) 

Before calculating the power spectra from our simulations, we add a component 
representing a point-mass inspiral matched to the beginning of our hydrodynamic 
merger wave form. This produces a spectrum with dE/df oc f 1 ^ 3 for point-mass 
inspiral at low frequencies. In essence, what we measure here is the number of 
orbits spent around a given frequency, weighted by the amplitude of the emission. 
We identify three frequencies of special interest. The frequency at which dynam- 
ical instability sets in is labeled fd yn - The frequency at the peak of the gravity 
wave luminosity is labeled f pea k- Last, the characteristic frequency of gravity wave 
emission for the merger remnant at late times in the calculation is labeled f osc . 

In Fig. 4, we show the gravity wave power spectrum for a Newtonian calculation 
with a synchronized initial condition, which includes radiation reaction effects but 
contains no 1PN terms. At low frequencies, we see the power-law behavior from the 
point-mass inspiral, with no contribution whatsoever from our calculated signal. At 
the dynamical instability frequency, we see a slight decrease in the gravity wave 
power, since the rapid plunge causes the binary to sweep up faster through a range 
of characteristic frequencies. The gravity wave power shows a plateau near the 
peak emission frequency, when the effects of the rapid infall are balanced by the 
large increase in gravity wave amplitude. At higher frequencies, there is another 
slight dip in emission, followed by a sharp peak marking the oscillation frequecncy 
of the merger remnant. With more careful handling of the late-time behavior of the 
system, we expect the peak to remain prominent, but our calculation most likely 
overemphasizes the coherence of the signal. 




FIGURE 4. Power spectrum calculated from a Newtonian coalescence calculation. The dotted 
and dashed lines represent the contributions of the point-mass inspiral and our calculated gravity 
wave signal for the hydrodynamic merger, respectively. The solid line represents the total power. 



A strikingly different power spectrum is obtained from our Post-Newtonian, ir- 
rotational merger, as shown in Fig. 5. The most significant difference between the 
two concerns the limit of dynamical instability. Newtonian gravity is strengthened 
by the addition of relativistic effects, which moves the dynamical stability limit 
to larger separations (and thus smaller frequencies). Additionally, the final plunge 
of the two NS is much faster, which significantly reduces the power in the region 
between fd yn and f pea k- At f pea k, the power is smaller than in the Newtonian case, 
but the signal is much more sharply defined. Similarly, the oscillation of the rem- 
nant at late times leaves a well defined imprint on the power spectrum, but the 
amplitude of the peak is approximately an order of magnitude lower than what is 
found in a Newtonian calculation. 

The frequencies of the two peaks seen in the spectrum, representing peak emission 
and the remnant oscillations, do give a strong clue to the nature of the NS EOS. 
While the frequency of peak oscillation is essentially the same in all our simulations, 
the width of the peak is seen to be strongly dependent on the EOS. The softer T = 2 
EOS shows a broad peak of emission in the frequency range / ~ 1500 — 3000 Hz, 
whereas the stiffer T = 3 EOS calculations have a peak much more focused around 
/ = 1800 — 2200 Hz, regardless of the initial spins. The remnant oscillations break 
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FIGURE 5. Power spectrum from a Post-Newtonian coalescence calculation. Conventions are 
as in the previous figure. 

this degeneracy. The oscillation frequency for the irrotational run is almost 15% 
less than that of the synchronized run with the same choice of EOS. The stiffer EOS 
also results in a more rapid oscillation than the softer one, although the frequencies 
are relatively similar. 

In general, several trends can be recognized regarding the strength of gravity wave 
emission from NS binary coalescence. The large dip at the dynamical instability 
limit is a general feature of PN calculations, regardless of the choice of initial spins 
or the EOS. It should be regarded as a consequence of the stronger gravity present 
in the PN systems. Additionally, PN simulations generally show similar amplitudes 
to their Newtonian counterparts near the characteristic frequency of peak emission, 
but significantly less power at higher frequencies, since the effect of strong gravity 
seems to be a quenching of gravity wave emission after the initial peak. Soft 
EOS generally show greatly reduced gravity wave emission at high frequencies, 
since they cannot support a stable, radiating, ellipsoidal configuration, and on 
relatively short timescales will produce a nearly spheroidal, non-radiating remnant. 
Finally, for binary systems with unequal-mass components, the magnitude of the 
gravity wave emission is strongly correlated with the mass ratio q [1,9]. Because 
the primary in such systems generally remains relatively undisturned, whereas the 
secondary is tidally disrupted and accreted onto the primary, a large component of 



the matter essentially does not contribute to the gravity wave signal. Thus, even 
if NS masses do not typically lie within a narrow range, there should be a strong 
bias observationally toward detection of nearly equal-mass systems. 
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